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ABSTRACT 

We study the hydrodynamical behavior occurring in the turbulent interaction 
zone of a fast moving red supergiant star, where the circumstellar and interstellar 
material collide. In this wind-interstellar medium collision, the familiar bow 
shock, contact discontinuity, and wind termination shock morphology forms, with 
localized instability development. 

Our model includes a detailed treatment of dust grains in the stellar wind, 
and takes into account the drag forces between dust and gas. The dust is treated 
as pressureless gas components binned per grainsize, for which we use ten rep- 
resentative grainsize bins. Our simulations allow to deduce how dust grains of 
varying sizes become distributed throughout the circumstellar medium. 

We show that smaller dust grains (radius < 0.045/im) tend to be strongly 
bound to the gas and therefore follow the gas density distribution closely, with 
intricate finestructure due to essentially hydrodynamical instabilities at the wind- 
related contact discontinuity. Larger grains which are more resistant to drag 
forces are shown to have their own unique dust distribution, with progressive 
deviations from the gas morphology. Specifically, small dust grains stay entirely 
within the zone bound by shocked wind material. The large grains are capable 
of leaving the shocked wind layer, and can penetrate into the shocked or even 
unshocked interstellar medium. Depending on how the number of dust grains 
varies with grainsize, this should leave a clear imprint in infrared observations of 
bowshocks of red supergiants and other evolved stars. 



Subject headings: Hydrodynamics — ISM: abundances — ISM: kinematics and 
dynamics — Stars: winds, outflows — Infrared: ISM 
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Introduction 



Hydro dynamical simulations of the circumstellar medium of low-to-intermediate mass 
stars (LIMS) have traditionally focused on the r norphology of the g as and neglected the 
presence of dust grains in the stellar wind (E.g. IWareing et al.l 120071). T hose simulations 
that do include dust tend to focus on the inner win d region (Woitkel 2006[) or deal w ith the 
proto- planetary disks surrounding very young stars (jPaardekooper fc Mellemal 120061 ). 



With recent advances in infrared observations, such as satellites like Spitzer and Her- 
schel, it has become possible to make detailed infrared observations of the circumstellar 
environment. The higher sensitivity of these satellites has enabled us to observe the weak 
signatures of the bowshocks that form when the wind of a moving star collides with the 
surrounding interstellar medium (ISM). Since dust grains are the primary source of infrared 
radiation, it becomes necessary to fully integrate them into the simulations to determine how 
the distribution of dust grains correlates with the morphology of the circumstellar gas and 
how the presence of dust influences the structure of the shock that occurs where the stellar 
wind interacts with the ISM. 

As an example, we have investigated the behavior of dust grains in the circumstellar 
environment of a fast moving red supergiant star. As input parameters we use the val- 
ues obtained for a-Orionis, which has been observed in infrared with the AKARI Satellite 



( lUeta et al.l 1200 81 ). This star shows evidence of a bow shock, indicating that it is moving 
through the ISM at a ve locity of iQn'J^'^'km/s with uh the dimensionless local hydrogen 
density JUeta et al.lboosh . 



So far, most simulations of fast-moving stars have focused on hot, high mass stars 
( Brighenti fc D'Ercole 1995 : Comeron fc Kaper 1998 : van Marie et al. 1 2006), wh i le sirn - 
ulations of the bowshocks of cool, lower mass stars were done by IWareing et al.l ( 120071 ). 
However, none of these simulations took the presence of dust into account. 



2. Modeling gas and dust interactions 

2.1. Governing equations 

If one wants to incorporate dustgrains in the hydro dynamical models, one should not 
only take into account the motion of the dustgrains themselves, it also becomes necessary 
to add the kinetic interaction between gas and dust to the equations of hydrodynamics. 



We use the MPI-AMRVAC code (iMehani et al, 



2007: 



Keppens et al 



20111 ). which solves 
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sets of near-conservation laws on an adaptive mesh. For this particular study, we added 
a new module to the code which calculates the behavior of dust grains as coupled to the 
familiar gas dynamic equations. The dust and gas are linked through drag forces which they 
exert on each other. Besides conservation of mass, the momentum and energy equations for 
the gas become 
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with p the gas density, v the gas velocity, p the thermal pressure and e the total energy 
density for the gas. We assume a standard ideal gas law for closure. The right hand side 
source terms for the energy evolution describe optically thin radiative cooling and the work 
done by the drag forces fd. Radiative losses depend on the hydrogen and electron particle 
densities (derived from p assuming full ionization with hydrogen mass mh) and involve a 
temperature dependent cooling curve A(T). As most of the gas in our simulation is at 
comparatively low temperature, we use a cooling curve A(T) that extends do wn to 1 K, with 



tabulated info based on numerical calculations done with the cloudy code ( iFerland et al. 



19981 ). The drag force for species d per unit volume fd is specified in Eq. |5l since we treat 
multiple dust species we sum over the individual drag forces. 



We follow the prescription from iPaardekooper fc Mellemal ( 120061 ). treating the dust as 
a gas without pressure, since the internal energy of a dustgrain only influences its surface 
temperature, but has no influence on its movement. Therefore, the motion of each dust 
species can be treated with the equations of conservation of mass and momentum: 

dpd 



dt 



^(PdVd) 

dt 



+ V • (pdVd) 
V ■ (pdVdVd) 
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fd 



(3) 
(4) 



where pd and Vd are the mass density and velocity of the dust species. The drag force fd is 
given by a combination of Epstein's d rag law for the subsonic regime and Stokes' drag law 
for the supersonic regime (lKwoklll975l ). which in closed form writes as 



— (1 — a)'Kndpa^A\-\/ A v 



(5) 



Its form involves the dust particle density rid, the particle radius ad (assuming that the dust 
particles are spherical), the velocity difference between dust and gas A v = v — Vd and the 
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thermal speed of the gas Vt = |a/3p/p (lKwoklll975l ). It also includes the sticking coefficient 
a, which is a measure of the percentage of atoms that stick to the dust grain after collision. 
This treatment neglects interactions between the different dust species. However, because 
of the low particle densities (<C lcm~^) and the small collision cross-section (< 1 /im^) of 
the individual grains, the g a s-dus t interaction will dominate over the interaction between 
dust particles. As in iKwokl (119751 ) we assume that a = 0.25, so 75% of the collisions are 
supposed to be elastic. 



2.2. Simulation parameters and setup 



For our simulation we use a 2D cylindrical grid with r = to 2 pc and z = —2 to 2pc. 
The base grid has a resolution of 80x160 gridpoints, and we allow grid adaptivity up to a 
maximum of six additional refinement levels, which translates in an effective resolution of 
5 120x10 240 gridpoints. The adaptive mesh is handled dynamically and traces the presence 
of dust, ensuring that the front of each expanding dust type is always fully resolved. We 
initialize the stellar wind by filling a spherical area of radius 0.1 pc around the origin with 
wind material. Since the star is moving through the ISM, we let the ISM flow past the star 
at a constant velocity. 

As input parameters for our simulations we use the observations of a -Qrionis, a typical 
example of a fast moving evolved star, as obtained by lUeta et al.l (j2008[ ). This gives us a 
gas mass-loss rate of Mgas = 3 x lO^^MQ/yr, injected at w ind velocity fwinri = 1 5 km/s. 
We estimate the dust mass-loss rate to be 2.5 x 10~^M^„, (|Verhoelst et al. Il2006[ l. Since 
the winds of evolved LIMS are believed to be dust-driven ( lKwoklll975l ). which means that 
the gas is dragged along with the grains, we assume that the dust and gas in the free- 
streaming wind are closely coupled, so that the velocity difference vanishes. Therefore, we 
initialize the dust grains with the same terminal velocity as the gas. This wind collides 
with an assumed homogeneous ISM, which has a particle density of nisM = 2 cm^^ and, in 
the frame of reference of the star, flows past the star with vi^m = 28.3 km/s. For the dust 
particles we assume a minimum radius of a^i^ = 0.005 /im and a maximum of ctmax = 0.25 fim 
( jPecin et al. II2006I ). with particle density distributed over the grainsizes as n{a) ~ a"^'^. 
We represent this distribution with particles of ten different sizes, taking n = 10 in Eq. 1-2, 
logarithmically distributed over the total grainsize interval. Each dust species effectively 
represents all dust grains within part of an interval A a^. The size of each grain type is given 
by 



log (fld) 



log (a„ 



(6) 
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We assume the dust to be com posed of silicates, for which the internal particle density is 
3.3 g/cm^ ( jPraine fc Lee Ill984j ). Together with the mentioned size distribution, this then 
fixes the dust particle densities rid, which in combination with the injection velocity (constant 
for all grain sizes) determine the dust mass-loss. The particle sizes and masses used in this 
paper can be found in Table [TJ 



3. Gas and dust distributions 

The results of our simulation are shown in Figs. [T]through|l]after a period of 75,000 years, 
which show the gas density and the particle densities for the individual dust species. 

In Fig. [Don the right we show the gas density, which clearly shows the bowshock and the 
wind termination shock, with the shocked gas layer in between. The distance b etween the 



shocks and the star is approximately 0.3 pc, which corresponds to observations ( lUeta et al. 



20081 ). The shocked gas region shows that the contact discontinuity between shocked wind 
and shocked ISM is subject to Rayleigh- Taylor type instabilities, which dominate the inter- 
action front. These instabilities result from the density difference between the shocked wind 
and ISM and are enhanced by the dust, which tend to move through the contact disconti- 
nuity. The temperature and absolute velocity of the gas are shown in Fig. [2J The shocked 
wind initially has a temperature of about 4,000 K. By the time it reaches the contact dis- 
continuity, the temperature has decreased to about 1 000 K due to radiative cooling. This 
decrease in temperature leads to compression as the thermal pressure of the gas decreases, 
resulting in a thin, high density feature at t he contact discontinuity . A m ore parametric 



study of this phenomenon has been shown in Ivan Marie fc Keppens I ( 120111 ). The shocked 



ISM has a higher temperature (~ 12, 000 K) because of the higher relative velocity of the 
ISM relative to the star (about a factor 2) and cools less efficiently due to its lower density. 

In the free-streaming wind, the dust grains have the same velocity as the gas, so they 
feel no drag force. As the dust grains cross the shock, they begin to experience a strong drag 
force because of the sudden difference in velocity with the gas. Because the dust initially 
has the same velocity as the pre-shock gas and the shock reduces the gas velocity by a 
factor 4 (according to the Rankine-Hugoniot conditions for an ideal gas), the initial velocity 
difference is | f„ind- The thermal velocity of the particles in the wind is of the same order of 
magnitude, so both the Stokes and Epstein drag forces (Eq. [5]) contribute to the total force. 
As the dust grains slow down due to the drag force, the Epstein regime starts to dominate. 
How quickly the velocity difference decreases depends on the size of the dust grains. Large 
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grains, which have a large momentum relative to their effective collision cross-section will 
take longer to slow down than small dust grains. As a result, the distribution of the dust 
grains over the gas becomes a function of the grain size. 

The smallest dust grains (left hand side of Fig. [1]) tend to follow the same pattern as the 
gas. As the wind slows down due to the wind termination shock, so do the smaller dust grains, 
which tend to pile up in the shocked wind region. For the smallest grains, the coupling to 
the gas is so strong that they follow the Rayleigh- Taylor instabilities (Fig. [T]). Some cross the 
contact discontinuity to enter the shocked ISM and are carried away downstream. This plot 
also allows us to determine which part of the gas is wind material and which is interstellar: 
the small dustgrains stay in the wind material. This is shown clearly in the region behind 
the star, where the contact discontinuity is almost parallel with the motion of the ISM. 
Figure 13] shows the distribution of larger dust grains (type 5, a5=0.0284/im), which have 
enough momentum to cross into the shocked ISM region and reach the forward shock. They 
are less susceptible to instabilities in the gas, though the effect of these instabilities is still 
visible. 

The largest dust grains can cross the shocked wind region due to their larger individual 
momentum and smaller collision surface relative to their mass. As a result, they mostly 
ignore the gas instabilities and eventually penetrate the forward bowshock to cross into the 
free-streaming (in the frame of the star) ISM. Since the free-streaming ISM is cold (assumed 
IK for this simulation), the Stokes' regime dominates over Epstein's drag law in Eq. |3 In 
the Stokes' regime, the drag force varies with the velocity difference squared. This makes 
the free-streaming ISM quite effective at stopping the oncoming dust grains (which move in 
the opposite direction, so the velocity difference is large) causing most of the particles to 
halt as soon as they cross into the free-streaming ISM (Fig. [3]). Only the largest particles 
can continue. These will penetrate deeply into the ISM as they are all but immune to the 
drag force. This is shown in Fig. HJ which shows the distribution of the type 8 dust grains 
(0.105/im) compared to the distribution of the smallest grains. The large grains almost 
completely ignore the morphology of the shocked gas shell and penetrate a considerable 
distance (~ 0.1 — 0.2 pc) into the unshocked ISM. Since they are very low in number, they 
will not influence the gas morphology. 

Since in nature the distribution of dust grains over the grain sizes is smooth, rather 
than existing of ten separate species, we can expect the distribution to be smoother than the 
relatively peaked positions found in our simulations. In the shocked gas, we could expect at 
least one, possibly several peaks in the grain density. The first peak would be in the shocked 
wind region, where the smallest grains will tend to pile up. A second peak may be visible 
just beyond the forward shock, where those grains that manage to cross the shocked gas 
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region encounter the unshocked ISM. In the unshocked ISM, the dust grains would spread 
out according to their size, with the largest grains extending furthest. This would not lead 
to any distinctive features. The tendency of large grains to pass through the shocked gas, 
while small grains are trapped, will deliver a disproportionately large number of large dust 
grains into the interstellar medium. For this particular simulation, the division lies by grains 
with a radius of approximately 0.045 /im. Smaller grains will stay inside the shocked gas, 
larger grains tend to continue outward; the distance reached by the individual dust species 
are quantified in Table [U 



4. Conclusions 

In the circumstellar nebula of a fast moving red supergiant, the small sized dust grains 
(a <0.045 /im) tend to follow the motion of the gas. Therefore, observations of infrared 
radiation emitted by these particles will be representative of the gas morphology as well as 
the morphology of the dust itself. The larger grains tend to create a morphology of their 
own, as they are only weakly coupled to the gas. 

Because the small dust grains will tend to dominate the observations as a result of 
their larger number and larger collective radiative surface area, the infrared observations 
of circumstellar nebulae can be used to determine the location of the shocked gas directly. 
However, observers should be aware of this phenomenon as it may influence observational 
results. For example, depending on the distribution of the number of grains over the grain 
sizes, a single shock may be observed as a series of separate features due to dust grains piling 
up in different dependent on their size. 

In the future we plan to extend our simulations to a more general parameter study and 
attempt to duplicate th e morphology of ob s erved circumstell ar nebulae around other cool 



stars, such as CW Leo (iLadjal et al. Il2010l : ICox et al. Il201ll ). We will also investigate in 



detail how the presence of dust influences the behavior of the gas. 
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Fig. 1. — This figure shows the gas density in g/cm^ (right) and particle density in ^/cm^ 
of the smallest dust grains (type 1, radius 0.005 fim) in the circumstellar medium of a 
fast-moving red supergiant star after 75,000 years of simulation time. The general shape of 
the bowshock between the stellar wind and the interstellar medium (ISM) has reached an 
equilibrium position. The interaction region shows strong Rayleigh- Taylor instabilities as a 
result of the density difference between the shocked wind and the shocked ISM. The dust 
tends to pile up at the contact discontinuity and follows the local instabilities. This also 
indicates which part of the gas is wind material as these small grains tend to pile up at the 
transition between wind and ISM. 



Fig. 2. — The gas velocity in km/s (left) and temperature in K (right) at the same time as 
Fig. [H The shocked ISM is hotter than the shocked wind, since it starts with more kinetic 
energy and cools less efficiently due to its lower density. Note that the velocity is not in the 
co-moving frame, but rather in the rest-frame of the ISM. 
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LoglO(number density) (dust n 5) 
-11 -10 -9 




Fig. 3. — Particle density in ^/cm^ of the intermediate sized dust grains (type 5, 
a=0.0284/im) from Table [H on the left compared to the gas density (right) at the same 
time as Figs. [T]and|2l These grains tend to follow the forward shock, rather than the contact 
discontinuity, but still show the influence of the instabilities in the gas. 
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0.5 pc 



Fig. 4. — The particle density distribution in #/cm^ of the larger dust grains (type 8, radius 
0.105 fim) compared with the smallest grains (type 1, radius 0.005 /xm), at the same time as 
Figs. [T]l3l These grains show little influence from the shocked gas layer and penetrate deeply 
into the unshocked ISM before being stopped by the drag force. 
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Table 1: Specification of tlie logaritlimic distribution of dust grain sizes. Column 2 gives 
the size of each particle species, column 3 the size interval which the species represent (see 
Eq.lH]), column 4 gives their mass and the last column the maximum distance that they reach 
in the direction of the motion of the star before being stopped by the drag force. 



name 


a[fim] 


A a [/im] 


mass per grain 


[g] 


max. distance [pc] 


1 


0.00500 


0.00136 


1.728x10" 


-18 


0.34 


2 


0.00772 


0.00346 


6.366x10- 


-18 


0.35 


3 


0.0119 


0.00535 


2.345x10- 


-17 


0.37 


4 


0.0184 


0.00826 


8.639x10- 


-17 


0.37 


5 


0.0284 


0.0128 


3.183x10- 


-16 


0.37 


6 


0.0439 


0.0197 


1.172x10- 


-15 


0.38 


7 


0.0679 


0.0304 


4.320x10" 


-15 


0.42 


8 


0.105 


0.0470 


1.591x10" 


-14 


0.53 


9 


0.162 


0.0726 


5.863x10- 


-14 


N/A'^ 


10 


0.250 


0.0441 


2.160x10" 


-13 


N/A'^ 



"This species has not been stopped by the drag force within the physical space of this simulation and continues 
to expand into the ISM 



